import scipy.stats
from scipy.interpolate import interp1d
import numpy
import math


#defined here to prevent circular import of spec.py and spectrum_maths.py
def _subtract(spec1, spec2):
    """
    Use spectrum_maths.subtract() instead.
    
    Returns a new spectrum whose counts are equal to that of spec1 - spec2.
    If spec2 has a different integration time to spec1, then it is scaled to 
    that of spec1 prior to subtraction.
    """
    
    if not isinstance(spec1, Spectrum):
            raise TypeError("Expected spectrum object for spec1 argument, got "+str(type(spec1))+" instead.")

    #only subtract a spectrum object 
    if not isinstance(spec2, Spectrum):
        raise TypeError("Cannot subtract object of type"+str(type(spec2))+" from Spectrum object")
    
    #only subtract spectra if they have the same wavelengths
    if not numpy.all(spec1.wavelengths == spec2.wavelengths):
        raise ValueError("Cannot subtract spectra with different wavelengths.")
    
    #only subtract spectra from the same spectrometer
    if spec1.spectrometer_id != spec2.spectrometer_id:
        raise ValueError("Cannot subtract spectra from different spectrometers"
                         "("+str(spec1.spectrometer_id)+' and ' + 
                         str(spec2.spectrometer_id) + ")")
    
    s =  Spectrum(spec1.filename, #keep spec1's filename
                  spec1.wavelengths, #wavelengths have not changed
                  spec1.counts - ((spec1.int_time / spec2.int_time)*spec2.counts),
                  spec1.capture_time,
                  spec1.spectrometer_id, #has not changed
                  spec1.int_time,
                  spec1.n_avg
                  )
    
    if spec1._is_dark is not None:
        s.set_is_dark(spec1._is_dark)
        
    return s


class Spectrum:
    def __init__(self, filename, wavelengths, counts, capture_time, spectrometer_id, 
                 int_time, n_avg, **kwargs):
        
        self.filename = filename #can be None
        self.wavelengths = numpy.array(wavelengths)
        self.counts = numpy.array(counts)
        self.capture_time = capture_time
        self.spectrometer_id = spectrometer_id
        self.int_time = int_time
        self.n_avg = n_avg
        self._is_dark = None
        
        for k,v in kwargs:
            if self.__dict__.has_key(k):
                raise ValueError("Cannot use a kwarg to override an existing attribute")
            self.__dict__[k] = v
 
    
    def __sub__(self, spec):
        """
        Subtract one spectrum from another using spec1 - spec2 syntax.
        """
        return _subtract(self, spec)
        
    def get_pixel_at(self, wavelength):
        """
        Returns a tuple of (pixel number, actual wavelength) for the pixel 
        closest to specified wavelength. If the specified wavelength is outside
        the range of the spectrum then raises ValueError.
        """
        if wavelength < self.wavelengths[0] or wavelength > self.wavelengths[-1]:
            raise ValueError("Requested wavelength is outside of spectrum range")
        
        idx = numpy.argmin(numpy.fabs(self.wavelengths - wavelength))
        return (idx, self.wavelengths[idx])
        
    
    def get_counts(self, wavelength):
        """
        Returns the counts at a specific wavelength. Linear interpolation will be used
        to retrieve counts at wavelengths in-between those that were recorded at.
        """
        if wavelength < self.wavelengths[0] or wavelength > self.wavelengths[-1]:
            raise ValueError("Requested wavelength is outside of spectrum range.")
        
        interpolator = interp1d(self.wavelengths, self.counts)
        
        return float(interpolator(wavelength))
        
    
    def set_is_dark(self, value):
        """
        Manually flag this spectrum as a dark spectrum (or not)
        """
        if not isinstance(value, bool):
            raise TypeError("set_is_dark expected boolean argument")
        self._is_dark = value
        
    
    def is_dark(self):
        
        #use the cached result if it is available
        if self._is_dark is not None:
            return self._is_dark
        
        #otherwise deduce whether this is a dark spectrum by comparing the CDF of the 
        #spectrum data to that of a uniform distribution using the Kolmogorov-Smirnov
        #test
        cumsum = numpy.cumsum(self.counts)
        cdf = cumsum/cumsum[-1]
        
        ks_result = scipy.stats.kstest(cdf, 'uniform')
        
        #compute the KS critical value at 99% confidence
        ks_critical_val = 1.628 / math.sqrt(cdf.shape[0])
        
        if ks_result[0] < ks_critical_val:
            self._is_dark = True
        else:
            self._is_dark = False
        
        return self._is_dark 
